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c/3 ! Solitary waves in a general nonlinear lattice are discussed, employing as a 

' model the nonlinear Schrodinger equation with a spatially periodic nonlinear 

p ^' coefficient. An asymptotic theory is developed for long solitary waves, that span 

O , a large number of lattice periods. In this limit, the allowed positions of solitary 

CZ3 I waves relative to the lattice, as well as their linear stability properties, hinge 

• ,-H ■ upon a certain recurrence relation which contains information beyond all orders 

^,^, of the usual two-scale perturbation expansion. It follows that only two such po- 

J^ I sitions are permissible, and of those two solitary waves, one is linearly stable and 

. *~^' the other unstable. For a cosine lattice, in particular, the two possible solitary 

waves are centered at a maximum or minimum of the lattice, with the former 
being stable, and the analytical predictions for the associated linear stability 
eigenvalues are in excellent agreement with numerical results. Furthermore, 
r^ . a countable set of multi-solitary-wave bound states are constructed analyti- 

cally. In spite of rather different physical settings, the exponential asymptotics 
approach followed here is strikingly similar to that taken in earlier studies of 



f"^ ' solitary wavepackets involving a periodic carrier and a slowly-varying envelope, 

which underscores the general value of this procedure for treating multi-scale 
solitary-wave problems. 
>■ 

•i-H 

X 

- - ■ 1 Introduction 

The nonlinear Schrodinger (NLS) equation is fundamental to wave propagation 
in fluid flows, optics, Bose~Einstein condensation (BEC) and applied mathe- 
matics [1, 2, 3, 4, 5, 6]. When a linear periodic potential (so-called linear lattice) 
is added, the resulting equation then models nonlinear beam transmission in 
an optical medium with a periodic transverse variation in the linear refractive 
index, and atom-atom interaction in Bose-Einstein condensates loaded in an 
optical lattice [7, 8, 9, 10]. Another interesting possibility, instead of a linear 
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lattice, is to allow the nonlinear coefficient of the NLS equation to vary peri- 
odically in space. In optics, such a periodic nonlinear coefficient (also called 
nonlinear lattice) would arise in the propagation of laser beams in a medium 
whose nonlinear refractive index is modulated in the transverse direction. The 
same model equation also applies to the study of Bose-Einstein condensates in 
a medium with a spatially-dependent scattering length. In optical and BEC 
experiments, a linear lattice can be created by beam interference [8, 11]. A 
nonlinear lattice in optics can be created by femtosecond-laser writing in fused 
silica [12], and in BEC such a lattice can be created by an optical Feshbach- 
resonance modulation of the scattering length [13, 14, 15]. Wave phenomena 
in linear lattices have been extensively studied in the past decade (see [10] for 
a review). In this paper, our interest centers on wave phenomena in nonlinear 
lattices. 

Wave phenomena in symmetric nonlinear lattices have been investigated in 
[16, 17, 18]. It was found that there exist solitary waves which are centered 
at a local maximum or minimum of the nonlinear lattice. Short solitary waves 
were found to be stable/unstable when centered at a local maximum/minimum 
of the lattice [17, 18], but the stability analysis for long solitary waves, that 
extend over a large number of lattice periods, was not conclusive [17]. Bound 
states comprising several fundamental solitary waves have also been reported 
numerically [16]. In addition to these studies in nonlinear lattices, work has 
also been done in the presence of both linear and nonlinear lattices [19, 20, 21] 
and in higher dimensions [22]. 

In this paper, we make an analytical study of long solitary waves and their 
linear stability properties in a general nonlinear lattice. Recognizing that the 
coupling between a long solitary wave and the relatively short nonlinear lattice 
is an exponentially small effect, beyond all orders of the usual multiple-scale per- 
turbation expansion in powers of the long-wave parameter, our analysis makes 
use of an exponential asymptotics technique. We show that long solitary waves 
can only be located at two positions in one period of the nonlinear lattice, 
regardless of the number of local maxima and minima in it. If the lattice is 
symmetric, the solitary-wave positions are simply the point of symmetry and 
half a lattice-period away from it; but for a general asymmetric lattice, these 
positions need to be determined by solving a certain recurrence relation that 
contains information beyond all orders of the long-wave parameter. Linear sta- 
bility analysis of long solitary waves is also performed: it follows from the same 
recurrence relation that one of the two solitary waves is linearly stable, while the 
other is unstable. For cosine lattices, in particular, the solitary wave centered 
at the maximum/minimum of the lattice is linearly stable/unstable, consistent 
with previous results [17, 18]. An analytical formula for the linear-stability 
eigenvalues is also derived; as expected, these eigenvalues are exponentially 
small in the long- wave parameter. The analytical predictions for the stability 
eigenvalues are compared with direct numerical results and excellent agreement 



is obtained. Lastly, we also show that an infinite number of multi-solitary- 
wave stationary bound states exist in the nonlinear lattice, and their analytical 
construction in terms of nonlocal solitary waves is presented. 

The exponential asymptotics procedure adopted in this paper closely resem- 
bles that used in recent studies [23, 24] of gap solitons in a linear lattice (see 
also [25, 26] for the first application of this technique to solitary wavepackets 
of the fifth-order KdV equation). The common thread in linear and nonlinear 
lattices is that the solitary wave is much longer than the period of the underly- 
ing lattice; hence, the coupling between these different scales is expected to be 
exponentially small, which invites similar exponential asymptotics treatment. 
As a result, the behavior of solitary waves in these rather different physical 
settings is remarkably similar. 



2 Preliminaries 

We study the NLS equation with a nonlinear lattice 

i^i + ^,, + (l+5(2/e))|^|'^ = 0, (2.1) 

where g{z/e) is a periodic function which describes the spatial variation of the 
nonlinear Kerr coefficient, and the parameter e > controls the length scale 
of this variation. Throughout this article, g{z/e) will be referred to as the 
nonlinear lattice. Eq. (2.1) is a model for the spatial propagation of a laser 
beam in a medium whose nonlinear refractive index is modulated periodically 
in the transverse direction (in this context, t is the direction of propagation), 
and for the dynamics of Bose-Einstein condensates whose scattering length (the 
counterpart of the nonlinear coefficient in (2.1)) changes periodically over space 
[4, 15]. Solitary-wave solutions in Eq. (2.1) and their linear-stability properties 
for even functions of g{z/e) were investigated in [17, 18]. In the long- wave 
limit (e ^ 1), profiles of solitary waves that span many lattice periods were 
determined by a multiscale perturbation analysis [17], but their stability was 
not ascertained. 

In this article, we analytically study long solitary waves (0 < e ^ 1) and 
their stability properties for general forms of the nonlinear lattice g{z/e). Specif- 
ically, denoting z/e = x, g{x) is assumed to be periodic with period d, 

g{x + d)= g{x) (2.2) 

for all real values of x. Also, without loss of generality, g{x) is taken to have 
zero X- average: 

i9) = lj g{x)dx = 0. (2.3) 



Solitary-wave solutions of Eq. (2.1) are sought in the form 

^{z,t) = ij{zy^'\ (2.4) 

where /U > is the propagation constant, and the real- valued function ip{z) 
solves the ordinary differential equation 

^zz + {l+g{z/e))ij^ -^liJ = ^, (2.5) 

under the boundary condition of V' — ^ as \z\ — )■ oo. 

Here, we will fix ^u = 0(1), and determine the solitary wave ip{z) as well as 
its linear stability for < e ^ 1. In this regime, the nonlinear lattice g{z/e) is 
rapidly varying and the wave profile ip{z), whose width is 0(1) for ^i = 0(1), 
spans many lattice sites. As e — )■ 0, g{z/e) features extremely rapid oscillations 
with zero mean and the effect of the nonlinear lattice on the solitary wave ip{z) 
can be dropped; thus, in this limit, ipiz) is expected to approach the familiar 
solitary-wave solution of the lattice-free NLS equation 

ip{z) ^ a sech. — - — , (2-6) 

where 

a = y^, p = 1/V^, (2.7) 

and zq denotes the location of the peak of the solitary wave 'ip{z). When 7^ e <C 
1, the nonlinear lattice g{z/e) will have a weak but non-negligible effect, and 
solitary waves bifurcate out from the limiting wave (2.6); this bifurcation will be 
the main focus of our investigation. It will be shown that such a bifurcation is 
possible only when zq takes two special values relative to the nonlinear lattice, 
resulting in two solitary- waves, out of which one is stable and the other unstable. 

It is worth mentioning that a different but equivalent analysis is to introduce 
scaled variables 

tp = eip, X = z/e, /I = e^fi, (2-8) 

so that Eq. (2.5) transforms into 

ipa^a^ + {l + g{x))i^^-l24^ = 0. (2.9) 

Thus, fixing fi = 0(1) and varying e <C 1 in our treatment above corresponds 
to fixing the lattice g{x) and varying /2 ^ 1 in Eq. (2.9). In this alternative 
treatment, the nonlinear lattice g{x) is no longer rapidly varying, but for ju <C 1, 
solitary waves ^(x) of Eq. (2.9) that bifurcate out from the edge /Iq = 
of the continuous spectrum (whose linear eigenmode is a constant) have low 
amplitude and are long compared to the lattice period. This type of solitary- 
wave bifurcation is akin to that studied in [23] with a linear lattice. In both 
treatments, however, the solitary wave is of long extent and spans many lattice 
sites, hence the fundamental feature of the problem remains the same. 



3 Multiscale perturbation solution 

We begin by reviewing the multiple-scale procedure followed in [17] for solving 
equation (2.5). The solution ip to this equation contains two scales, reflected 
in the 'slow' variable z and 'fast' variable x = z/e. Introducing explicitly these 
variables by writing ip = ^(x,z), Eq. (2.5) then becomes 

dl + ^d^d, + ^d^] V + (1 + g{x))i;^ -fii; = 0. (3.1) 

Now we expand 'ip{x,z) into the two-scale perturbation series, 

ip{x, z) = i^oix, z) + eV'i(x, z) + e^ilj2{x, z) -\ . (3.2) 

Substituting this expansion into Eq. (3.1) and from various orders of e, we get 
the following hierarchy of equations, 

-dli^o = , (3.3a) 

-dl^Pi = 2a^9,Vo , (3.3b) 

-^l^P2 = dli,^ -h 29^5,Vi + (1 + 9{x))^l - /i^o , (3.3c) 

-92^3 = ^'V'l + 1d:,d,i,2 + 3(1 + 5(2;))V'oV'i - /^V'l , (3.3d) 

-dli,^ = dl^2 + 2d^d,ij^ + 3(1 + g{x))il,li^2 + 3(1 + g{x))i^oi^l - fiiJ2 . 

(3.3e) 

From Eq. (3.3a) and the requirement that ipo be bounded, we get 

i^o = Mz)- (3.4) 

Substituting this equation into (3.3b) and requiring ipi to be bounded, we get 

V'i = V^i(z). (3.5) 

Equations (3.3c)-(3.3e) are forced linear equations and can be written in the 
unified form 

-dllpn = Qn{x,z), (3.6) 

where Qn depends on ipj with j < n. The x-dependence of Qn derives from the 
nonlinear lattice g{x) which is d-periodic. Thus Qn is d-periodic in x and can 
be expanded into a Fourier series, 



00 



Qn{x,z)= V g„(z)e2-'"^/'^. (3.7) 



m=— 00 



We require solutions ipn to be also d-periodic in x. The necessary and sufficient 
condition for the existence of such a solution is that the constant Fourier mode 
qo in (3.7) vanish, 

qo = {Qn)=0. (3.8) 



Here (•) is the average with respect to x as defined in (2.3). Then the solution 
ipn to Eq. (3.6) is 

1pn = -d~^Qnix,z)+'lpniz) , (3.9) 



where 

d^'Qn=Er-^) '/™(^)e^~/^ (3.10) 
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and Ipn is a function of z which is determined by the solvabihty condition (3.8) 
for the equation governing V'n+2 (see below). 

Next we determine tpoiz) and tpiiz). Specifically, to obtain ipo{z), we return 
to equation (3.3c). Substituting (3.5) into (3.3c) and recalling that g{x) has 
zero average (see (2.3)), the solvability condition of this equation gives 

92vio + ^il-/i^o = 0, (3.11) 

whose solution is 

■00 = ^(z) = asech — - — , (3-12) 

where a and /3 are defined in (2.7), and zq is a constant. This result agrees with 
(2.6), as anticipated earlier. 

To determine the solution ipi{z), we consider equation (3.3d). Using the 
above expression for tpQ and the zero average of g{x), the solvability condition 
for (3.3d) gives 

{d^/dz^ -fi + 3A'^)7jji = 0, (3.13) 

whose solution is 

V^i = C^'(^), (3-14) 

where (^ is a constant. It is clear that ipi simply shifts the location of the center 
of the leading-order term, ipQ = A{z), in the perturbation expansion (3.2). To 
remove the ambiguity in the position of ipQ, we require that V'l be orthogonal 
to ijjQ, hence C = 0, and 

Mz) = 0. (3.15) 

Now we determine the solution iIj2 to Eq. (3.3c). Utilizing solutions f/'o smd 
ipi in Eqs. (3.12) and (3.15), '(/'2 can be written as 

Mx, z) = -[d-^g{x)]A'{z) + Mz) . (3.16) 

When this solution is inserted into equation (3.3e), the solvability condition of 
this equation gives 

{d^/dz^ - H + 3A^) i)2 = -3a A^ , (3.17) 

where 

« ^ {{d^^g?) > 0. (3.18) 



The bounded solution i{j2 to this equation is 

Mz) = a{A^ - 2a'^A). (3.19) 

Combining (3.4), (3.5), (3.15), (3.16) and (3.19), the perturbation series 
solution for the solitary wave ip{x, z) of Eq. (2.5) takes the form 

V'(x,z) = A{z)+e^\-[d^^g{x)\ A^{z)+a[A^{z) -2a^A{z)] \+0{e^), (3.20) 



where A{z) is given by (3.12). 

Notice that the location zq of the peak of the function A{z) in Eq. (3.20) is 
arbitrary at this stage, since equation (3.11) which determines A{z) is transla- 
tion invariant. However, the original equation (2.5) for the solitary wave is not 
translation invariant due to the presence of the nonlinear lattice g'(x), and it is 
unlikely that the solitary wave can be arbitrarily located relative to this nonlin- 
ear lattice. Indeed, a very similar situation arises in a linear lattice [6, 23, 27], 
and there it was shown that solitary waves can only be located at two positions 
relative to the lattice. In the present problem, the result turns out to be similar: 
only two values of zq are permissible for truly localized solitary waves in a non- 
linear lattice, as will be established in the next section by utilizing exponential 
asymptotics. 

Before ending this section, it should be pointed out that truly localized 
solitary waves ip{z) of Eq. (2.5) must satisfy 

g'{z/e)i)^{z)dz = 0. (3.21) 

This constraint can be obtained by multiplying Eq. (2.5) with V''(-2) and then 
integrating from — oo to oo. An analogous constraint was noted in [27] for the 
linear lattice problem, and it was pointed out that, if the lattice is symmetric, 
this constraint predicts only two possible locations for truly localized solitary 
waves — one at the point of symmetry and the other half a period away from 
it [6, 27]. But for general asymmetric linear lattices, it does not seem feasible 
to determine the locations of solitary waves based on this constraint alone [23] . 
Similarly, in the problem at hand, if the nonlinear lattice g{z/e) is symmetric, 
the constraint (3.21) can also predict the two locations of true solitary waves; if 
the lattice is asymmetric, however, this approach would again fail. The reason 
is that, when the perturbation series solution (3.20) is substituted into the 
constraint (3.21), all terms in the series make contributions of the same order 
of magnitude to the integral in (3.21). Hence it does not seem possible to solve 
for zq without having obtained all terms in the perturbation series (3.20). The 
same difficulty also appears in the linear lattice problem [6, 23] and suggests 
the need for a perturbation theory beyond all orders. This task is taken up 
below by employing an exponential-asymptotics procedure in the wavenumber 
domain [23, 25]. 



4 Growing tails of exponentially small amplitude 

In this section, we determine the location of the sohtary wave ip{z) by the 
exponential asymptotics method. Our approach closely resembles that for gap 
solitons in a linear lattice [23], thus only the key ideas and steps will be given. 
For further details, we refer the reader to [23] (see also [25]). 

The ensuing analysis is based on the fact that, if ip{z) given by (3.20) is 
required to decay upstream (z <^ —1)) then this solution of Eq. (2.5) will 
contain a growing tail downstream (z ^ 1) for generic values of the position zq 
of the solitary wave core. Specifically, if the upstream asymptotics of ip{z), as 
given by the leading-order term in the perturbation expansion (3.20), is 

V'~2ae(^-^o)/'^, z^-oo, (4.1) 

then the downstream asymptotics of the solution will be 

V'~2ae-(^"^")/'^ + /7e(^~^")/'^, 2>1, (4.2) 

where H{e, zq) is the growing-tail amplitude. As we will show later, H vanishes 
only at two special values of zq (relative to the nonlinear lattice), thus only 
two truly-localized solitary waves exist. The tail amplitude H{e,zo) turns out 
to be exponentially small in e, so this growing tail can not be captured by the 
perturbation series (3.20) and has to be obtained by carrying this expansion in 
powers of e beyond all orders. 

Following the exponential asymptotics procedure in the wavenumber domain 
[23], we introduce the Fourier transform of ip{x,z) with respect to the slow 
variable z, 

^ 1 r°^ 

ip{x,K) = — ^{x,z)e-'^'dz. (4.3) 
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Substituting the perturbation series solution (3.20) for iIj{x, z) into this Fourier 
transform, we find that 



,i^) = ^e-^^Osech(|/?K)|l 



2„2 



-{l+P^K^)u{x)-2a 



+ 



(4.4) 

where i^{x) = a — d^'^g{x). This expansion in the wavenumber domain is 
disordered when K = 0(l/e), suggesting that V has pole singularities a,t K = 
0(l/e). The residues of these singularities are exponentially small due to the 
exponentially small value of the sech function in (4.4) aX K = 0(l/e). As 
we will show later, these singularities of exponentially-small residue contribute 
exponentially-small but growing tails in the physical solution '4){x, z). Thus the 
main goal is to determine the locations and residues of pole singularities in ^. 



To this end, we replace the disordered expansion (4.4) by a uniformly-valid 
expression, 

^ = e-*-^^Osech(|/3i^jC7(x,K), (4.5) 

where eK = k. We then take the Fourier transform of Eq. (3.1) with respect 

to z, 

i)^^ + 2iK^^^ - K^i^ + £^(1 + 5(x))V^ - eV ^ = , (4.6) 

and upon substituting Eq. (4.5) into (4.6), we obtain the following equation for 
U 

U^x + 2iKUx - K^U - e^nU 

+ (l+g(x))cosh—K dX ^ ' — ^ dp ^ ',, ,V ^ =0- 

^ ^ ^^ 2 J_^ cosh^i^ 7-00 '^cosh^i^ cosh 2M 

(4.7) 



4.1 The recurrence equation 

In the limit e — )■ and for k away from singularities of U{x,k.), the main 
contribution to the double integral in Eq. (4.7) comes from < A < «; when 
K > and from k < A < when k < 0. The integral equation (4.7) then 
simplifies to [23, 25] 



Ua;x + 2iKU^-K'^U + 4(l + g(x)) I dXU(x,K-X) I dpU(x, X- p)U(x, p) = . 



-g{x)) I dXU{x,K-X) [ dpU{x,X-p)Uix,p) 
Jo Jo 

(4.8) 

The solution to this simplified integral equation can be posed as a power series 
in K, 

U{x,K) = ^Y.Un{x)K''. (4.9) 

n=0 

Substituting this power series into (4.8), we obtain the following recurrence 
equation for Un, 

, o =Un- 2l^ 

dx^ dx 

- 2(1 + g{x)) V Vn-m ,~X[\ Y] UrU^.rrl{m - r)! . (4.10) 

m=0 ^ ' r=0 

To be consistent with (4.4), we set 

[/o = 1, Ux= 0, (4.11) 

which serve as the initial conditions for the recurrence iteration (4.10). Note 
that this recurrence system does not involve p and e; it only depends on the 
functional form of the periodic lattice g(x). 



4.2 Behavior near singularities 

When K is near the singularities of U{x, k), the reduced integral equation (4.8) 
does not apply. We now examine the behavior of U{x, k) near its singularities, 
based on the original integral equation (4.7). 

These singularities occur at k « kq, where the linear part of Eq. (4.7) is 
satisfied, i.e., 

U^^ + 2moU^-KlU = 0. (4.12) 

The bounded solution to this linear equation is C/ ~ g-««;ox Since U is expected 
to be d-periodic in x in view of (4.5), therefore, kq = ziz2ir/d. Singularities near 
kq = zizin/d, ziz6ir/d, . . . are also possible, but they are subdoniinant and will 
not be considered. To avoid ambiguity, we denote 

Ko = 2TT/d, (4.13) 

and examine singularities at k ~ ±«;o- 

In order to determine the behavior of the solution U near the singularity 
K ~ kq, we first introduce the 'inner' wavenumber 

{ = ^i^, (4.14) 



that is, K = kq + eS, with ^ = 0(1). Guided by the analysis in [23, 25], we 
expand 

„—iKox 

U = ^^{<^o{0 + e^f{x,0 + ---). (4.15) 

The reason for the leading-order term in this expansion being 0(e~^) will be 
explained later (see Eq. (4.22)). 

Now we examine the integral equation (4.7) near the singularity k ~ kq. The 
dominant contributions to the double integral in (4.7) come from the following 
regions: (i) A ~ and p ~ 0, (ii) A ~ k with p ^ and with p ^ k. Taking 
into account the leading-order term near k ~ kq in (4-15) and the leading order 
term near k ~ in (4.9), we can calculate these dominant contributions, and 
Eq. (4.7) near k ~ kq yields 

(4.16) 
where U = e~*''°^[/. Substituting the expansion (4.15) into (4.16), the terms 
of 0(e~^) and 0(e~^) are automatically balanced. At 0(e~^) we have 

dlf = fMO + f^MO - 3(1 + 9{x)) r oje^^'^/' csch ^M^ -u;)doj, 

(4.17) 
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and the solvability condition requires that the x-average of the right-hand side 
of Eq. (4.17) equals to zero: 

(m + f)MO - 3 r oje^^'^/^ csch ^$o(e -u;)du; = 0. (4.18) 

Recalling that /? = 1/ yjjl, this integral equation is identical to that encountered 
earlier in our analysis of gap solitons in a linear lattice [23]. Hence, its solution 

is [23, 25] 

^oiO = -r^TT^ I -^cjy{s)e-'P^ds, (4.19) 



1 + /32^2 _/^± sin^s 
where 



^/ 2 cos^s 3scoss\ 
(«) = C( — T + — - - — ^ ) ' (4-20) 



\sins sms sin s 

the contours C extend from to iticx) for Im(^) < and Im(^) > respectively, 
and C is a complex constant which will be determined later. Note that the 
function <5o(^) given by (4.19) is analytic everywhere in the complex plane C, 
save for two simple poles at ^ = ibi//3. Also, <&o('?) satisfies the integral equation 
(4.18) only in the complex ^-plane outside the strip — 1//3 < Im(^) < 1//3, as 
was explained in [23]. 

The behavior of the function ^o(0 near its singularities ^ = ±i//3 is 

*«tt)--Tff^ tt-^i). (4.21) 

and its large-^ asymptotics is [23, 25] 

^o(0^^|4 (e^^)- (4-22) 

This fourth-order algebraic decay rate of <l>o('?) a-t large ^ implies the order e~^ 
in Eq. (4.15). From Eq. (4.21), we then obtain the behavior of ip{x,K) near 
K = Ko/e =F V/3 as 



/93/7 „-iKQ(x+xo) ■ 



From the symmetry of the Fourier transform for real functions, ip is also singular 
at i^ = — Ko/e =F i//^, and 
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4.3 The growing tail and locations of solitary waves 

We now take the inverse Fourier transform of ip{x, K) in order to calculate the 
growing tails in tlj{x, z) due to the pole singularities (4.23)- (4.24). The inverse 
Fourier transform is 

il,{x,z)= I i){x,K)e'^''dK, (4.25) 

where the contour C is taken along the line '^{K) = — 1//3 and to pass below 
the poles ai K = ±Ko/e — V/^ (the reason for this choice of the contour is 
explained in [23]). Also, the contour C passes above the pole of sech(7r/3i^/2) at 
K = — i//3, to be consistent with the desired upstream behavior of tp in (4.1). 

Indeed, for z ^ — 1, the dominant contribution to the inverse Fourier trans- 
form (4.25) comes from the pole of sech(7r/3i^/2) at -ftT = — i//3 and one obtains 
the upstream solution (4.1). On the other hand, for z ^ 1, both the singulari- 
ties (4.23)~ (4.24) si K = ±Ko/e - i//3 and the singularity of sech(7r/3E:/2) at 
K = i/ (3 contribute, and one obtains the downstream solution behavior as 



(.-.„)//3 ^ 47r^^_^^,„/2, ^.^^^^^^ _ ^^^(,_,„)/;3 



V' ~ 2ae-^^-^"^/'^ + ^^ Q-^P^ol^^ sin(KoXo - 0)e^^-^o^/'^ (^ > 1) , (4.26) 

where C > and 6 are the amplitude and phase of the complex constant C, 

C = Ce'^ . (4.27) 

Equation (4.26) is one of the key results in this paper. It shows that a growing 
tail of exponentially small amplitude appears far downstream in the solution 
tp{x,z), except when 

sin(KoXo -9) = 0, (4.28) 



I.e., 



xo = e/Ko , (6 + tt)/ko , (4.29) 



Thus, exactly two solitary waves, located at these values of xq, are obtained in 
the nonlinear lattice equation (2.5). 

4.4 Determination of the constant C 

It remains to determine the complex constant C. This constant cannot be 
obtained from the local analysis around the singularities k ~ ±kq, but it can be 
computed by solving the recurrence relation (4.10). For this purpose, we derive 
the asymptotics of the recurrence functions [/„ for large n, which depends on C. 
The derivation is based on the idea that the 'inner' solution (4.15) of U{x,k) 
near the singularities, when (^ — t- oo, must match the 'outer' solution (4.9) of 
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U{x,k) away from the singularities. First, from the inner expansion (4.15) and 
the asymptotics (4.22), we see that 

-\2C 1 

f/(x,K)~^- ^6-*"°^ (e<|K-Ko|<l), (4.30) 

5 (k - Ko) 



or 



12C 1 

5 (k - kqY 



U ^ -J [cos(ko2; — 9) — isiTi.{nQX — 9)] (e ^ |k — kq] <C 1). 



(4.31) 

From the symmetry relation U{x,k.) = U*{x,—k*) for real functions ^(x,z), 
we also have 

12C 1 

U ~ —^7 Y4 [cos(Koa; - 6*) + isin(Koa; - 9)] (e ^ |k + kqI ^ !)• 

5 (k + kq) 

(4.32) 
These two asymptotic expressions can be combined and re-written as 



1924c 

U ^-' 



5(ac2 - 4) 



2\4 



(e< |k±ko| <. 1) 



cos(Koa^ — 0) — i — sin(Koa: — 6 

(4.33) 

This expression is consistent with (4.31)-(4.32) near the singularities. More 
importantly, it has the property that, when expanded into power series of k, 
the coefficients of all even powers of k are purely real and the coefficients of all 
odd powers of k are purely imaginary, as required by the outer solution (4.9), 
with Un given by the recurrence relation (4.10) under the initial conditions 
(4.11). 

Expanding (4.33) into power series of k. and requiring this expansion to be 
consistent with the outer solution (4.9), we obtain the asymptotic behavior of 
Un for n » 1 as 

TTl Tfl 

U2m ~ ^^;^ COs(koX - e) , U2m+1 ~ "« ^ 2m+l sin^KQX - 9) , (4.34) 

where the coefficient D is related to C by 

a = ^K>. (4.35) 

Thus, by solving the recurrence relation (4.10) and from its large-n asymptotics, 
we can obtain D and 9, and hence the complex constant C. Since the recurrence 
equation (4.10) depends only on the lattice function g{x), the constant C also 
depends only on the lattice g{x) and not on ^, e. 
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5 Linear stability problem 

In this section, we determine the hnear stabihty of the two sohtary waves whose 
locations are given by Eq. (4.29). We will show that the solitary wave located at 
xo = / kq is linearly stable, while that located at {O + i:)/ kq is linearly unstable, 
and the unstable eigenvalue is exponentially small in e. This calculation follows 
the approach used in [23, 26] for the stability analysis of solitary wavepackets 
of the fifth-order KdV equation and gap solitons of the NLS equation with a 
linear lattice. 

Let V's(-2) be a solitary wave of Eq. (2.5), whose leading-order term -00 = 
A{z) in (3.12) is centered at zq = zqs-, where zqs = exos, and xqs is one of the 
two positions given in Eq. (4.29). Consider the perturbed solution 



■^izd) 



''^* Up,{z) + [v{z) + tf (z)]e^* + [v*{z) - w*{z)]e^**\ , (5.1) 



where v,w <^ 1 are normal-mode perturbations. Substituting this perturbed 
solution into Eq. (2.1) and neglecting nonlinear terms in {v,w), we obtain the 
linear-stability eigenvalue problem [6] 

LqLiv = -X^v. (5.2) 

where 

d^ d^ 

Lo = ^ + {I + g{x))i^l - fi, Li = ^ + 3{1 + g{x))il^l - fi. (5.3) 

If there exist eigenvalues A with positive real parts, then the solitary wave is 
linearly unstable. Otherwise it is linearly stable. 

When e ^ 1, the discrete eigenvalue A is small. To calculate this eigenvalue, 
we expand v into powers of A, 

v = vo + X'^vi + X'^V2 + ■■■ . (5.4) 

Substituting this expansion into Eq. (5.2), at 0(1), we have 

LoLiVo = 0. (5.5) 

Recall that 

Loi^{x,z;xo) = 0, (5.6) 

where ip{x, z; xq) is the nonlocal solution to Eq. (2.5) whose leading term ipQ = 
A{z) in (3.12) is centered at z = exQ. Taking the derivative of this equation 
with respect to xq and then setting xq = xqs, we get 



L — 
dxo 



= 0. (5.7) 

Xo=XOs 
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Hence the solution to Eq. (5.5) is 



^0 « 



dtp 
vo = 

Notice from Eq. (4.26) that 
^„-(.-.o.)//3 , ^ 



(5.8) 



y^ ^ '^^-i^-^os)/^ + 4^^o^^g-./3.o/2. cos(koXo. - e)e(^--o^)/^ z » 1, 



(5.9) 

thus vq contains a growing tail and is nonlocal. Since the discrete eigenfunction 
V must be localized, this growing tail in vq must be cancelled by another growing 
tail in X'^vi. 

The equation for vi is found from Eq. (5.2) at O(A^) as 

LqLiVi = -Vq. (5.10) 

Letting Livi = wq, we first solve 

LqWq = -Vo- (5.11) 

Since the growing tail in vq is exponentially small, its contribution to a likewise 
exponentially-small growing tail in vi (through Eq. (5.10)) can be ignored, 
because the localized (algebraically small) terms in vq will turn out to create 
a relatively larger (i.e., algebraically small) growing tail in vi. Thus, in the 
calculation of wq, it suffices to take vq as (5.8), with V' given by the pertur- 
bation series (3.20). The corresponding solution wq can be expanded into a 
perturbation series 

wo = eB{z) + e^wo{x,z) + --- . (5.12) 

Inserting this expansion into (5.11), at 0(e) we obtain 

-dl wo = B"{z) - fiB{z) + (1 + g{x))A\z)B{z) - A'{z) . (5.13) 

The solvability condition of this equation yields the governing equation for B(z) 



as 

i2 



d'B „ .2 



dz2 

hence 



fiB + A^B = A'{z), (5.14) 



B{z) = ^iz-zos)Aiz). (5.15) 

Now we solve for vi from LiVi = wq. This solution can be expanded as 

vi=eF{z) + e^vi{x,z) + ... . (5.16) 

Substituting this expansion into LiVi = wq, at 0(e) we get 

-d^vi = F"{z) - fiF{z) + 3(1 + g{x))A\z)F{z) - ^(z - zos)A{z) . (5.17) 
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The solvability condition of this equation gives 

^-liF + 3A^F = ^{z-zos)A. (5.18) 

Note that the homogeneous solution A'{z) of (5.18) is not orthogonal to the 
inhomogeneous term, hence the solution F(z) to Eq. (5.18) is nonlocal. For 
our purpose, we seek a solution F{z) so that F{z) — )■ as z — )■ — oo and 
F{z) —7- FQe^^~^'-'"'''^ for z S> 1, where Fq is a constant. To determine Fq, we 
multiply Eq. (5.18) by A'(z) and then integrate from — oo to z S> 1, 



A'iz 



F"{2) - nF{z) + 3A^{z)F{2) 



di = - / {z-zos)A{z)A'{z)dz. 



(5.19) 
Employing integration by parts on the left-hand side and using the large-z 
asymptotics of F{z) above, we can obtain the left integral in terms of Fq. The 
right integral approaches a constant as z — )■ +oo, which can be readily obtained 
using the functional form (3.12) of A{z). After these calculations, Eq. (5.19) 
yields Fq = afS'^/S. Hence, the corresponding large-z asymptotics of vi is 

v^ ~ ^ee("-""=)/^ (z > 1) . (5.20) 

8 

Inserting this growing tail of vi and the growing tail of vq in (5.9) into the 
expansion (5.4) of v and utilizing the relation (4.35), the eigenvalue formula is 
then found to be 

5 g-7r/3Ko/2e 

A2 = --K^7r/3D -g cos(Koa;os - ^) , (5.21) 

where D and 6 are obtained by solving the recurrence equation (4.10). Note 
that this eigenvalue is exponentially small in e. In addition, since D > 0, the 
solitary wave located at xqs = ^//^o is linearly stable, and the one located at 
xqs = {0 + tt)/kq is linearly unstable. 

We remark in passing that the eigenvalue A given by formula (5.21) is two 
orders larger (in e) than the eigenvalues found in earlier studies for solitary 
wavepackets of the fifth-order KdV equation [26] and Bloch-wavepacket solitons 
of the NLS equation with a linear lattice [23]. However, this does not imply 
that the linear instability (for one of the two solitons) in the present nonlinear 
lattice problem is stronger than those in the fifth-order KdV equation and the 
NLS equation with a linear lattice. The reason is that the small parameter e 
in [23, 26] measures the wave peak amplitude, while in the present analysis e 
is a long-wave parameter. Indeed, the peak amplitude of solitary waves in our 
analysis is 0(1) rather than 0(e) (see Eq. (2.6)). In the end of Sec. 2, we 
mentioned an alternative treatment where, through rescaling, the solitary wave 
becomes long and also features low amplitude of 0(e). This equivalent analysis 
is the proper counterpart of those in [23, 26]. After the rescaling (2.8), in fact, 
the stability eigenvalue is of the same order in e as in [23, 26]. 
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6 Numerical results 

In this section, we present numerical results for solitary-wave profiles and their 
linear-stability eigenvalues, and make a comparison with the above analytical 
results. The numerical algorithms for these computations can be found in [6]. 
In our computations, we take /i = 1 and the nonlinear periodic lattice to be 

9{x)=gocosx, (6.1) 

with go = 1. In this case, d = 2tt and kq = 2TT/d = 1. The numerical 
procedure for solving the recurrence relation (4.10) is similar to that in [23]. 
Our computation confirms that Un for n 3> 1 indeed approaches the asymptotic 
form (4.34) with 

B = 0.1638, 6' = 0. (6.2) 

Thus our theory predicts that the two solitary waves are located at xq = 
(maximum of g{x)) and xq = vr (minimum of g{x)), or equivalently zq = and 
zq = eir respectively. To borrow the terminology of gap solitons in linear lattices, 
we call the solitary wave located at the lattice-maximum zq = 'on-site', and 
the other one located at the lattice- minimum zq = evr 'off-site'. Numerically, 
we have computed these two solitary waves (for each value of e), and found 
them to be indeed located at the two zq positions. To demonstrate, these 
solitary waves for e = 0.5 and 0.15 are displayed in Fig. 1. Notice that the 
on-site solitary waves have a single hump, while the off-site ones have double 
humps. In addition, when e is small, both on-site and off-site solitary waves 
have approximately a sech profile, in agreement with the perturbation series 
solution (3.20). 

Next we numerically determine the linear stability of these solitary waves 
and make a comparison with our analytical results. The whole linear-stability 
spectra for the on-site and off-site solitary waves in Fig. l(b,d) at e = 0.15 
are shown in Fig. 2(a,b). The spectrum in Fig. 2(a) lies entirely on the 
imaginary axis, indicating that this on-site solitary wave is linearly stable. In 
this spectrum, a pair of purely imaginary discrete eigenvalues can be seen. 
These are the counterparts of our analytical imaginary eigenvalues given by 
Eq. (5.21) with xqs = O/kq = 0. The spectrum in Fig. 2(b) contains a 
real positive eigenvalue, indicating that the underlying off-site solitary wave is 
linearly unstable, in agreement with our analytical prediction. In particular, 
this positive eigenvalue is the counterpart of our analytical positive eigenvalue 
given by Eq. (5.21) with xqs = {0 + '^)/no = vr. 

Now we quantitatively compare the numerical linear-stability eigenvalues 
with the analytical formula (5.21). For this purpose, we have numerically ob- 
tained the discrete eigenvalues (as those in Fig. 2(a,b)) for on-site and off-site 
solitary waves at various values of e, and the results are shown in Fig. 2(c,d) 
by dotted lines. The analytical eigenvalues from the formulae (5.21) for the 
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Figure 1: Profiles of on-site and off-site solitary waves in the nonlinear lattice 
equation (2.5) at two values of e = 0.5 and 0.15 (with fi = 1). Upper row: on- 
site solutions; lower row: off-site solutions; left column: e = 0.5; right column: 
e = 0.15. The vertical stripes mark locations of high nonlinear-lattice values. 

on-site and off-site cases are also plotted as solid lines, and excellent agreement 
can be seen for both cases. This verifies that the analytical formula (5.21) is 
asymptotically accurate. 

Finally, we numerically examine how these solitary waves evolve nonlinearly 
under weak perturbations. For this purpose, we again consider the on-site and 
off-site solitary waves in Fig. l(b,d) at e = 0.15. We perturb these waves 
initially by a small phase gradient as 



^(z,0) = V(2)e' 



I'JZ 



(6.3) 



where ip{z) is the solitary wave and 7 = 0.01. This phase-gradient perturbation 
gives the solitary wave a small initial 'push'. Evolutions of the on-site and off- 
site solitary waves under this perturbation are obtained by simulating Eq. (2.1) 
with the above initial condition (6.3), and the results are displayed in Fig. 3. 
It is seen that the on-site solitary wave is not affected by this perturbation and 
stays at its initial on-site position (see Fig. 3(a)), consistent with the linear 
stability of this on-site solitary wave established in Fig. 2(a). On the other 
hand, under the same perturbation, the off-site solitary wave moves from its 
initial off-site position to a nearby on-site position and then oscillates around 
it (see Fig. 3(b)). When t — )■ 00, the oscillation eventually dies out, and the 
solution evolves into a stationary on-site solitary wave. This nonlinear evolution 
scenario is consistent with the linear instability of this off-site solitary wave, as 
shown in Fig. 2(b). 
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Figure 2: (a,b) Stability spectra for the on-site and off-site solitary waves in 
Fig. l(b,d) respectively (e = 0.15); (c,d) comparison between numerical and 
analytical discrete eigenvalues of on-site and off-site solitary waves at various 
values of e; dotted lines: numerical values; solid lines: analytical values from 
formulae (5.21). The arrows mark locations of the value of e in (a,b). 





Figure 3: Nonlinear evolutions of the on-site (a) and off-site (b) solitary waves 
in Fig. l(b,d) under phase-gradient perturbations (6.3). Red bars represent 
locations of high nonlinear-lattice values. 

7 Bound states 

In addition to the two fundamental solitary waves obtained earlier, the nonlinear- 
lattice equation (2.5) also admits higher-order solitary waves, or bound states 
[16]. These can be analytically constructed by matching the tails of more than 
one of the nonlocal solitary waves discussed in Sec. 4. A similar construction 
has been detailed in [24] for bound states of the NLS equation with a linear 
sinusoidal periodic potential (see also [25]). Here we shall sketch the analysis 
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for bound states involving two nonlocal solitary waves in a nonlinear lattice. 

For simplicity, we assume the symmetric cosine nonlinear lattice (6.1), with 
given depth gQ. In this case, kq = 1 and 9 = 0. Then we consider two nonlocal 
solitary waves of (see Sec. 4), ip^{z) and ^^(z), whose main 'sech' humps 
are centered at Zq = ex^ and z^ = —exQ, respectively, with Zq > 0, and 
Zq + Zq ^ 1 (large separation). In addition, let ip^{z) — t- as z — )■ ±00, so 
that the right-hand tail of ip^iz) and the left-hand tail of tp'^{z) are nonlocal. 
According to (4.26), the right-hand tail of ip~{z) is given by 

For symmetric lattice functions g{x), Eq. (2.5) is invariant with respect to 
reflection in z (z — t- —z). In addition, if ip{z) is a solution to (2.5), so is —ijj{z). 
Thus the left-hand tail of ip^{z) can be obtained from (4.26) after reflection in 

z as 

V^+(z) ~ ±2ae(^-^tt)//3 ^ ^^l^^-^N2e gi^^+ g-(.-4)//3 (^ « 4) . (7.2) 

Here, the upper sign in (7.2) corresponds to the case where the main humps 
of the two nonlocal waves have the same polarity (sign), while the lower sign 
pertains to the case of opposite polarity. To obtain a solitary wave (bound 
state) comprising these two nonlocal waves, we require that the right-hand tail 
(7.1) of '4^~{z) and the left-hand tail (7.2) of ^'"''(•2) match smoothly in the 
overlap region —Zq <C z <C Z(^. This requirement gives 



e^ 



sm Xq = sm Xq 



27rC/3^ 



,7r/3/2e g~e(x++x„ )//3 _ (7^3^ 



These matching equations are identical to those in [24] after a scaling in 
{xq ,Xq ,e,C), and their solutions can be taken directly from [24]. Specifically, 
for fixed /U > 0, e > and given sign (polarity), these equations admit an 
infinite number of solutions {xq,x^). Each solution corresponds to a bound 
state whose leading-order approximation is 

'(/'(z) ~ asech — —^ it a sech — w^j (7-4) 

and a continuous family of bound states is obtained when e or n varies. Each 
bound-state family contains triple branches and, for fixed e, these branches dis- 
appear when n falls below a certain threshold (or equivalently, for fixed /x, these 
branches disappear when e falls below a certain threshold). To demonstrate, 
we take the cosine nonlinear potential (6.1) with qq = 1 and e = 0.15. In this 
nonlinear lattice, a family of bound states comprising two fundamental solitons 
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of the same polarity is numerically obtained and displayed in Fig. 4. The power 
curve -P(/u) of this family, defined as 



^(/^) 



7p'^{z;n)dz, 



(7.5) 



contains three branches. On the lower branch, the bound state comprises two 
on-site fundamental solitons which are separated approximately by 8 lattice 
sites. On the upper branch, the bound state comprises two off-site fundamental 
solitons which are separated approximately by 7 lattice sites. On the middle 
branch, the bound state comprises an on-site and an off-site fundamental soliton 
which are separated approximately by 7.5 lattice sites. These solution branches 
exist only when n > fic ^ 0.7887; thus, they do not bifurcate from infinitesimal 
linear waves at the edge ^u = of the continuous spectrum of Eq. (2.1). These 
numerical results are in good agreement with the theoretical predictions based 
on the formulae (7.3)-(7.4). 






Figure 4: A family of bound states with same polarity in Eq. (2.5) for the 
cosine nonlinear potential (6.1) with go = ^ and e = 0.15: (a) the power curve; 
(b,c,d) bound-state profiles at points of the power branches marked by the same 
letters in (a). 



8 Concluding remarks 

In this article, we developed an asymptotic theory for long solitary waves and 
their linear stability in a general nonlinear lattice. Based on exponential asymp- 
totics, we showed that long solitary waves can only be located at two positions 
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relative to the nonlinear lattice, regardless of the number of local maxima and 
minima in the lattice. In general, these positions are determined by a certain 
recurrence relation that includes information beyond all orders of the usual 
multiple-scale perturbation expansion. From the same recurrence relation, one 
may also deduce that, of these two solitary waves, one is linearly stable and the 
other is unstable. If the lattice is symmetric, then the solitary-wave positions 
are simply the point of symmetry and half a lattice-period away from it. In 
particular, for the special cosine lattice, the solitary wave centered at the max- 
imum/minimum of the lattice is linearly stable/unstable. We also derived an 
analytical formula for the linear-stability eigenvalues, which are exponentially 
small with respect to the long-wave parameter (the ratio between the lattice 
period and the width of the solitary wave). The predictions of this analytical 
formula were compared against numerical results and excellent agreement was 
observed. Finally, it was pointed out that an infinite number of multi-solitary- 
wave bound states are possible in a nonlinear lattice, and their analytical con- 
struction was presented. 

The exponential asymptotics procedure used in this investigation closely 
resembles that in [23, 24] for linear lattices and in [25, 26] for the fifth-order 
KdV equation. In fact, the integral equation (4.18), which plays a key role in 
the analysis, actually arises in all these three different physical models. This 
suggests that our asymptotic procedure in the wavenumber domain is a possibly 
universal treatment of multiscale solitary-wave problems, and it is likely to find 
applications in other physical settings as well. 
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